M-SSA

Рассмотрим временной ряд динамики средней цены (доллар/фунт) на грейпфруты в США в Цетрально-Западном и Северо-Восточном регионе за двадцатитрехтилетний период (длина ряда 276 точек) и его периодограмму.

data.mw <- read.csv('APU0200711411.csv', sep = ',', header = T)
data.ne <- read.csv('APU0100711411.csv', sep = ',', header = T)
df <- data.frame(data.mw$APU0200711411[1:276], data.ne$APU0100711411[1:276])
colnames(df) <- c("Midwest", "Northeast")
price <- ts(df, start = c(1980, 1), end = c(2002,12), frequency = 12)
plot(price)

spec.pgram(price, log='no', detrend = FALSE, taper = 0)

Графики и периодограммы очень схожи. Неудивительно, ведь это цена на одни и те же фрукты в разных частях страны. Так как сигнал условно одинаковый, дополнительная информация из второго канала (регион Northeast) может помочь лучше отделиться от шума. Сравним разложение цены Midwest методом SSA и M-SSA.

s.mw <- ssa(price[,1])
plot(wcor(s.mw, groups = 1:20))

plot(s.mw, type = "vectors", idx = 1:15)

plot(s.mw, type = "paired", idx = 1:15)

Видим, что 2-3 компонента перемешалась с трендом.

Попытаемся улучшить разделимость внутри сигнальных компонент с помощью DerivSSA.

fs.mw <- fossa(s.mw, nested.groups = c(1:7))
plot(wcor(fs.mw, groups = 1:20))

plot(fs.mw, type = "vectors", idx = 1:15)

plot(fs.mw, type = "paired", idx = 1:15)

Разделимость стала намного лучше.

Посмотрим на восстановленный сигнал и acf остатков.

r.mw <- reconstruct(fs.mw, groups = list(Trend = c(3:7), 
                                        Seasonality = c(1:2)))
plot(r.mw, add.residuals = FALSE)

acf(residuals(r.mw))

По acf остатков можем предположить, что это шум.

Теперь применим Multi-channel SSA.

s <- ssa(price, kind = "mssa")
plot(wcor(s, groups = 1:20))

plot(s, type = "vectors", idx = 1:15)

plot(s, type = "paired", idx = 1:15)

Компоненты, в том числе полугодовая гармоника, хорошо отделились от шума, тренд не перемешан с гармониками сигнала.

Посмотрим на восстановление и сравним полученные разложения.

r <- reconstruct(s, groups = list(Trend = c(1, 4, 5, 8, 9),
                                  Seasonality = c(2:3, 6:7)))
plot(r, add.residuals = FALSE,
     plot.method = "xyplot",
     superpose = TRUE, auto.key = list(columns = 3))

plot(r$Trend[,1], col = "green", ylab = "Trend (Midwest)")
lines(r.mw$Trend, col = "red")
legend("topleft", 
       legend = c("MSSA", "DerivSSA"),
       col = c("green", "red"), lty = 1)

plot(r$Seasonality[,1], col = "green", ylab = "Seasonality (Midwest)", ylim = c(-0.15, 0.17))
lines(r.mw$Seasonality, col = "red")
legend("topleft", 
       legend = c("MSSA", "DerivSSA"),
       col = c("green", "red"), lty = 1)

spec.mw <- spec.pgram(residuals(r.mw), log='no', detrend = FALSE, plot = FALSE)
spec <- spec.pgram(residuals(r)[,1], log='no', detrend = FALSE, plot = FALSE)
plot(residuals(r.mw), col = "red",  ylab = "Residuals (Midwest)", ylim = c(-0.15, 0.17))
lines(residuals(r)[,1], col = "green")
legend("topleft", 
       legend = c("MSSA", "DerivSSA"),
       col = c("green", "red"), lty = 1)

plot(spec.mw$freq, spec.mw$spec, col = "red", type = "l",
     xlab = "frequency", ylab = "spectrum", main = "Периодограмма остатков")
lines(spec$freq, spec$spec, col = "green",)
legend("topright", 
       legend = c("MSSA", "DerivSSA"),
       col = c("green", "red"), lty = 1)

Результаты схожи, однако применение M-SSA позволило нам отделить компоненты полугодовой сезонности от шума.

2D-SSA

Загрузим изображение, добавим к нему текстуру и посмотрим, сможет ли 2D-SSA выделить изображение.

raw <- readRaw('image.raw')
matrix <- matrix(as.numeric(raw$fileRaw), nrow = 400, ncol = 400, byrow = TRUE)
x <- 1:(400*400)
noise <- matrix(5*sin(x)+10*cos(2*x), nrow = 400, ncol = 400, byrow = TRUE)
image <- matrix + noise
s <- ssa(image, kind = "2d-ssa", L = c(25, 25))
plot(s, type = "vectors", idx = 1:20,
     cuts = 255, layout = c(10, 2),
     plot.contrib = FALSE)

plot(wcor(s, groups = 1:30), scales = list(at = c(10, 20, 30)))

r <- reconstruct(s, groups = list(Noise = c(7:8, 14:15)))
plot(r, cuts = 255, layout = c(3, 1))

Искусственно наложенный шум хорошо отделился.